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Abstract. Higher moments of event-by-event net-proton multiplicity distributions are applied 
to search for the QCD critical point in the heavy ion collisions. It has been demonstrated that 
higher moments as well as moment products are sensitive to the correlation length and directly 
connected to the thermodynamic susceptibilities computed in the Lattice QCD and Hadron 
Resonance Gas (HRG) model. In this paper, we will present measurements for kurtosis (k), 
skewness (S) and variance (a 2 ) of net-proton multiplicity distributions at the mid- rapidity 
(\y\ < 0.5) and 0.4 < p T < 0.8 GeV/c for Au+Au collisions at N /s^"=19.6, 39, 62.4, 130 and 
200 GeV, Cu+Cu collisions at v /s7^7=22.4, 62.4 and 200 GeV, d+Au collisions at v /s7^7=200 
GeV and p+p collisions at ^/s NN =62.4 and 200 GeV. The moment products Ka 2 and Sa of net- 
proton distributions, which are related to volume independent baryon number susceptibility 
ratio, are compared to the Lattice QCD and HRG model calculations. The Ka 2 and Sa of 
net-proton distributions are consistent with Lattice QCD and HRG model calculations at high 
energy, which support the thermalization of the colliding system. Deviations of Ka 2 and Sa for 
the Au+Au collisions at low energies from HRG model calculations are also observed. 



1. Introduction 

One of the main goals of heavy ion collision is to explore the phase structure of hot, dense nuclear 
matter [1]. Finite temperature Lattice QCD calculations demonstrate that with vanishing 
baryon chemical potential (hb = 0), the phase transition from the hadronic phase to the Quark 
Gluon Plasma (QGP) phase is a smooth crossover [2]. The corresponding transition temperature 
is about 170 — 190 MeV [3111]. At large jib region, the phase transition between hadronic phase 
and QGP phase is of first order [5] and with a second order end point at the boundary towards 
the crossover region, which is so called QCD Critical Point (CP) [6]. Although many efforts 
have been made by theorist and experimentalist to search for the CP, its location or even the 
existence is still not confirmed yet [7j- The first principle Lattice QCD calculation at finite ^b 
are difficult due to the fermion sign problem. Several techniques, such as Re-weighting, Image 
baryon chemical potential and Taylor expansion [8] etc., have been developed to overcome those 
problems and make the Lattice QCD calculable at finite \ib region. However, large uncertainties 
are still there. 

Experimentally, heavy ion collision provides us a good opportunity to search for the CP. To 
access the QCD phase diagram, we can tune the chemical freeze out temperature (T) and baryon 
chemical potential (hb) by varying the colliding energy. As the characteristic signatures of the 



CP in a static and infinite medium is the divergence of the correlation length (£) and increase 
of non-Gaussian fluctuations. Thus, non-monotonic signals of CP are expected to be observed 
if the evolution trajectory (T,^lb ) in the QCD phase diagram of the system pass nearby the 
critical region and the signals are not washed out by the expansion of the colliding system. 

Due to the finite size effect, rapid expansion and critical slowing down etc., the typical 
correlation length (£) developed in the heavy ion collision near the QCD critical point is a small 
value about 2 — 3 f m [9]. Recently, model calculations reveal that higher moments of conserved 
quantities distributions are proportional to the higher power of the correlation length [W\ 
such as fourth order cumulant < (<5iV) 4 > -3 < (5N) 2 > 2 ~ £ 7 , where 6N = N — M, N is the 
particle multiplicity in one event and M is the averaged particle multiplicity of the event sample. 
On the other hand, the higher moments as well as moment products of conserved quantities 
distributions are also directly connected to the corresponding susceptibilities in Lattice QCD 
[12j [13] and HRG model [2] calculations, for e.g. the third order susceptibility of baryon 
number (xe?) is related to the third cumulant (< (5Nb) 3 >) of baryon number distributions 
as x^b = < O^Yb) 3 >/VT 3 ; V,T are volume and temperature of system respectively. It has 
been found that the volume independence baryon number susceptibility ratio can be directly 
connected to the moment products of baryon number distributions as kg 1 = XbVxb ano - 

(3) (2) 

So = Xb IXb ' w hich allows us to compare the theoretical calculations with experimental 
measurements. Theoretical calculations also demonstrate that the experimental measurable 
net-proton number (proton minus anti-proton number) event-by-event fluctuations can reflect 
the baryon number and charge fluctuations |15j . Thus, higher moments of the net-proton 
multiplicity distributions are applied to search for the QCD critical point in the heavy ion 
collisions [16} \T7\ 118] . When approaching the QCD critical point, the moment products kg 1 
and Sa will show large deviation from its Poisson statistical value. The skewness is expected to 
change its sign when system evolution trajectory in the phase diagram cross phase boundary |19j . 

In year 2010, RHIC Bean Energy Scan (BES) program [20] was carried out to map the first 
order QCD phase boundary and search for the QCD critical point by tuning the colliding energy 
from 39 GeV down to 7.7 GeV with the corresponding fig coverage about 100 — 410 MeV. With 
the large uniform acceptance and good capability of particle identification STAR detector, it 
provides us very good opportunities to find the QCD critical point with sensitive observable, if 
the existence of QCD critical point is true. 

2. Experimental Method 

The data presented in the paper are obtained using the Solenoidal Tracker at RHIC (STAR). 
Those are Au+Au collisions at ^/s^=19.6 (year 2001), 39, 62.4 (year 2004), 130 and 200 
GeV (year 2004), Cu+Cu collisions at ^/s^=22.4, 62.4, 200 GeV, d+Au at y^=200 GeV 
(year 2003) and p+p collisions at y/s^=62A (year 2006), 200 GeV (year 2009). The main 
subsystem used in this analysis is a large, uniform acceptance cylindrical Time Projection 
Chamber (TPC) covering a pseudo-rapidity range of \rj\ < 1 and 2n azimuthal coverage. As a 
primary tracking device, it can measure the trajectories and momenta of particles with transverse 
momenta above 0.15 GeV/c. To ensure the purity and similar efficiency, the protons and anti- 
protons are identified with the ionization energy loss {dE/dx) measured by TPC of STAR 
detector within 0.4 < pt < 0.8 GeV/c and mid-rapidity (\y\ < 0.5). Several track quality cuts 
are also used to select the tracks with good quality in each event. We require the distance 
of closest approach (dca) to the primary vertex of proton (antiproton) track less than 1 cm 
to suppress the contamination from secondary protons (antiproton). To study the centrality 
dependence of higher moments, centralities are determined by the uncorrected charged particle 
multiplicities (dN c h/dr]) within pseudo-rapidity \rj\ < 0.5 measured by TPC. By comparing 
measured dN c h/dr] with the Monto Carlo Glauber model results, we can obtain the average 



number of participant N part , impact parameter b and number of binary nucleon-nucleon collisions 
Nun for each centrality class. 



3. Analysis Method 

In this section, we will introduce you how to perform the higher moments analysis with 
experimental data. First, we will show you the definition of the cumulants and various moments 
used in our analysis, such as standard deviation (cr), skewness (S) and kurtosis (k). Then, we 
will discuss about the Centrality Bin Width Effect (CB WE). 



3.1. Moments and Cumulants of Event-by-Event Fluctuations 

In statistics, probability distribution functions can be characterized by the various moments, 
such as mean (M), variance (a 2 ), skewness (S) and kurtosis (k). Before introducing the above 
moments used in our analysis, we would like to define cumulants, which are alternative methods 
to the moments of a distribution. The cumulants determine the moments in the sense that any 
two probability distributions whose cumulants are identical will have identical moments as well. 

Experimentally, we measure net-proton number event-by-event wise, N p ^p = N p — Np, which 
is proton number minus antiproton number, in the mid-rapidity (|y| < 0.5) and within the 
transverse momentum 0.4 < px < 0.8 GeV/c. In the following, we use N to represent the 
net-proton number N p -p in one event. The average value over whole event ensemble is denoted 
by < N >, where the single angle brackets are used to indicate ensemble average of an event- 
by-event distributions. 

The deviation of N from its mean value are defined by 

SN = N- < N > . (1) 
Then, we can define the various order cumulants of event-by-event distributions as: 

Ci, N =< N >, C 2 ,n =< (SN) 2 >, C 3:N =< (5N) 3 >, C A , N =< [8N) A > -3 < (5N) 2 > 2 . (2) 

An important property of the cumulants is their additivity for independent variables. If X 
and Y are two independent random variables, then we have Cix+Y = Cj x + Cj y for ith order 
cumulant. This property will be used in our study. 

Once we have the definition of cumulants, various moments can be denoted as: 

\C2,N) 1 {C2,N) 

And also, the moments product ko 2 and So can be expressed in term of cumulant ratio. 



2 C± N C 3yN 
ko = - — ,So = - — . (4) 

With above definition of various moments, we can calculate various moments and moment 
products with the measured event-by-event net-proton fluctuations for each centrality. 



3.2. Centrality Bin Width Effect (CBWE) Correction 

The centralities in this analysis are determined by the uncorrected charged particle multiplicity 
(N c h) measured at middle pseudo-rapidity (\r]\ < 0.5) by the TPC of the STAR detector. Before 
calculating various moments of net-proton distributions for one centrality, such as 0—5%, 5—10% 
, we should consider the so called Centrality Bin Width Effect (CBWE) arising from the impact 



parameter fluctuations due to the finite centrality bin width. This effect must be corrected, as 
it may cause different centrality dependence. To formulate and demonstrate the centrality bin 
width effect, we write the event-by-event net-proton distributions in one centrality: 

P(N) = J2"if {i) (N),C£u> i = l), (5) 

i i 

where the u>i and (N) are the weighted and net-proton distributions for ith impact parameter 
in one centrality, respectively. From eqs. (2)-(5) and (11), we can calculate the various order 
cumulants for the distribution P(N) as below: 

Ci,jv = Y,u H Cl N = J2"i<N> i (6) 



C 2 ,N = (X) ^ C 2,n) + C 2,C\ N ( 7 ) 
i 

C 3 , N = ( W ^3,7v) + C z,c\ N + 3 x C i,q iJV ,i,q iJV ( 8 ) 

i 

C4,n = (y2uiCl N )+C 4C i +4xC 1C i 1C i 

i 

~\~ 6 X C-, ( r ii \2 1 f^ ^ [C-i t^ti )(Ci r*i 1 ) 

1 '( Lj 1,N> ' i '°2,iV 1 '°1,JV i '°l,JV' 1 '°2,iV 

where the C\ N (k = 1,2,3,4) are the k th order cumulant for net-proton distribution f^(N); 
the C' k x (X = C % m N , k,m = 1, 2, 3, 4) are k th order cutmulant for random variable X = N 
under the discrete probability distribution Prob{X)=uji] the C 1X1Y =< XY > — < X >< Y > 
(X = C l k N , k = 1, 2, 3, 4) are the first order joint cumulant for random variable X, Y under the 
discrete probability distribution Prob(X, Y)=u)i. We find that the higher order cumulants C^^n 
(k = 1,2,3,4) can be expressed by the addition of two parts, one is the weighted average of 
the same order cumulant of each sub-distribution /W(7V), and the other part is the cumulant of 
lower order cumulant under the discrete weighted distributions, which stems from the fluctuation 
of impact parameters within the centrality and results in the CBWE. 

Experimentally, the smallest centrality bin is determined by a single uncorrected reference 
multiplicity value measured by TPC. Generally, we usually report our results for a wider 
centrality bin, such as — 5%, 5 — 10%, ...etc., to supress the statistical fluctuations. To eliminate 
the centrality bin width effect, we calculate the various moment for each single N c h within one 
wider centrality bin and weighted averaged by the number of events in each N^. 
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K = — =y^0J r K r , (12) 



where the n r is the number of events in r th refmult and the corresponding weight u r = n r /Y2 n r . 
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4. Results 

In this section, we will present beam energy and system size dependence for the various 
moments (M,a, S, k) as well as moment products (ko~ 2 ,So~) of net-proton distributions for 
Au+Au collisions at ^/s^ = 19.6, 39, 62.4, 130, 200 GeV, Cu+Cu collisions at ^/s^ = 22.4, 
62.4, 200 GeV, d+Au collisions at ^/s^ = 200 GeV and p+p collisions at ^/s^ = 62.4, 200 GeV. 
First, we will show the typical event-by-event net-proton multiplicity distributions from different 
colliding systems. Then we studied the centrality as well as energy dependence of various 
moments and moment products. The systematic errors are estimated by varying the following 
requirements for p(p) tracks: DCA, track quality reflected by the number of ot points used in 
track reconstruction and the dE/dx selection criteria for p(p) identification. The statistical and 
systematic error are shown separately by lines and brackets, respectively. 

4-1- Event-by-Event Net-proton Multiplicity Distributions 

Event-by-event net-proton multiplicity distributions for Au+Au collisions at ^/s NN = 39 GeV 
measured within 0.4 < px < 0.8 GeV/c and \y\ < 0.5 are shown in Fig. [TJ Going from peripheral 
to central collisions, it is found that the distributions become wider and more symmetric for 
central collisions. 
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Figure 1. Typical event-by-event net-proton multiplicity distributions for Au+Au collisions at 



/s^ = 39 GeV. 



4-2. Centrality Dependence of Moments and Moment Products of Net-proton Distributions 
The centrality (N part ) dependence for various moments (M, a, S, k) of net-proton multiplicity 
distributions from Au+Au collisions at ^/s NN =39, 62.4 and 200 GeV, Cu+Cu collisions at 
^/s NN =22.4, 62.4 and 200 GeV, d+Au collisions at ^/s NN =200 GeV and p+p collisions at 
^/s NN =62.4 and 200 GeV are shown in Fig. [2] and Fig. [3j We find that M and a increase with 
N part monotonically, while S and n decrease with N part . The centrality and energy dependence 
of S and k indicate that the net-proton distributions become more symmetric for more central 
collision and higher energies. 

The dashed lines in the Fig. [2] and Fig. [3] represent the expectations from Central Limit 
Theorem (CLT) when assuming the superposition of many identical and independent particle 
emission sources in the system [TTl [18] . In Fig. [2] and US the centrality dependence of various 
moments can be well described by the dashed lines expected from CLT. Especially in Fig. [3l 
the various moments of p+p and d+Au collisions follow the CLT lines of Cu+Cu collision at the 
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Figure 2. Centrality dependence of various moments of Figure 3. Centrality dependence of various moments of 

net-proton multiplicity distributions for Au+Au collisions net-proton multiplicity distributions for Cu+Cu collisions 

at ^/s^ =39, 62.4, 200 GeV. The dashed lines shown at ,/s 
in the figure are expectation lines from Central Limit 
Theorem (CLT). 



NNI = 22.4, 62.4 and 200 GeV, d+Au collisions at 
y/s NN = 200 GeV and p+p collisions at ^/s NN = 62.4 
and 200 GeV. The dashed lines shown in the figure are 
expectation lines from Central Limit Theorem (CLT). 



corresponding energy very well. This also supports the identical independent emission sources 
assumption. Fig. [3] shows the centrality dependence of moment products Sa and na 2 of net- 
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Figure 4. Moment products (na 2 and Sa) of net-proton distributions for Au+Au, Cu+Cu, 
d+Au and p+p collisions. 



proton distributions, which are directly related to the baryon number susceptibility ratio in 

Lattice QCD and HRG models as kg 2 = XbVxb an d Sa = XbVxb ' f° r P+P' Cu+Cu and 
Au+Au collisions at various colliding energies. Sa shows a weak increase with centrality, while 
the na 2 shows no centrality dependence. 



4-3. Energy Dependence of the Moment Products (Sa and Ka 2 ) 
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Figure 5. Energy dependence of moment products (na 2 
and Sa) of net-proton distributions for central Au+Au 
collisions (0 - 5%, 19.6 GeV: - 10%, 130 GeV: - 6%). 
The red dashed lines denote the HRG model calculations, 
in which the Sa = tanh(/Us/T) and na 2 = l. The 
empty markers denote the results calculated from Lattice 
QCD @]. 



Figure 6. Energy dependence of moment products (na 2 
and Sa) of net-proton distributions for Cu+Cu central 
collisions (0 - 10%, 22.4 GeV: - 5%). The dashed lines 
shown in the figures are from HRG model calculations, 
in which Sa = tanh(^/r) and na 2 =l. The fis and T 
values are taken from |21| . 



In Fig. [5l we show the energy dependence of the Sa and na 2 of net-proton distributions 
for most central Au+Au collisions (0 - 5%, 19.6 GeV: - 10%, 130 GeV: - 6%). Lattice 
QCD [4| and HRG model [13] calculations are also shown for comparison. Lattice QCD results 
are obtained with time extent N T = 6 and phase transition temperature at /Ub=0, T c = 175 
MeV. The red dashed lines of the HRG model in the upper panel and lower panels are evaluated 
by Sa = tanh(//B /T) and na 2 = 1, respectively, where the hb/T ratio at chemical freeze-out 
is parameterized as a function of colliding energy based on reference [22] . The corresponding 
baryon chemical potential {^b) at chemical freeze-out for each energy is shown in the upper 
band of the Fig. [5j We find that the moment products (na 2 and Sa) of Au+Au collisions at 
- v /s NN =200, 130, 62.4 GeV are consistent with Lattice QCD and HRG model calculations. The 
Sa and Ka 2 for Au+Au collisions at ^s NN = 39 GeV deviates from HRG model calculations. 
Surprisingly, na 2 from Lattice QCD calculations at ^/s NN = 19.6 GeV show a negative value [3J. 
However, due to the limited statistics, the statistical errors of experimental data are large at 19.6 
GeV. The STAR is running to take more data at 19.6 GeV in year 2011. Those deviations could 
be linked to the chiral phase transition [23] and presence of QCD critical point [24J. Recent 
linear a model calculations demonstrate that the forth order cumulant of the fluctuations for 
a field will be universally negative, when the QCD critical point is approached from cross-over 
side [23]. It will cause the measured na 2 as well as kurtosis (k) of net-proton distributions to 
be smaller than their Poisson expectation values. 

Fig. [6] shows the energy dependence of na 2 and Sa for Cu+Cu central collisions. The 
red dashed lines in the figure are obtained from the HRG model by using the formula 
Sa = tanh()UB /T) , where the hb and T are from thermal model fits of the particle ratios. 
We find that our experimental data is consistent with HRG model expectations for Sa of net- 
proton distributions. While the Ka 2 deviates from HRG model calculations and monotonically 



decrease as the collision energy decreases. 



4-4- Charged Particle Density (dN c h/dr]) Scaling of Sa and Evidence of Thermalization in 
Heavy Ion Collisions 
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Figure 7. Sa of net-proton distributions as a function 
of charged particle density at mid-rapidity (dN c h/dn) for 
various colliding systems. The dashed lines in the figures 
is the fitting lines. 



l - 



eo 10 -. 



10 - 



Lattice QCD HRG 

O N^.T -175 MeV 



jA 



STAR Preliminary 



STAR Data 

A d+Au 200 GeV T P+P 200 GeV 

Cu+Cu 22.4 GeV * Au+Au 62.4 GeV 

O Cu+Cu 62.4 GeV A Au+Au 130 GeV 

0-5% •Cu+Cu 200 GeV * Au+Au 200 GeV 



10" 



10" 



u /T 

B 



10 



Figure 8. Sa of net-proton distributions as a function of 
baryon chemical potential over temperature ratio (hb/T) 
for most central collisions (0—5%, 130 GeV: 0—6%, d+Au 
: — 20%). The red dashed line in the figures is the result 
of the HRG model. Lattice QCD results with iV T = 6 and 
T c = 175 MeV are also shown in the figure. 



The Sa of net-proton distributions for various colliding systems including Au+Au collisions 
at ^/J^ = 62.4 and 200 GeV, Cu+Cu collisions at ^/s^ = 22.4, 62.4 and 200 GeV, d+Au 
and p+p collisions at ^/s NN = 200 GeV, as a function of dN c h/drj are shown in the Fig. [7] 
with double logarithm axis. It is obvious that for a fixed colliding energy, such as ^/s NN =200 
GeV, the moment products of Sa of net-proton distributions for different system size the p+p, 
d+Au, Cu+Cu and Au+Au collisions have power law dependence on the charged particle density 
{dNch/dn). Thus, we fit the Sa of net-proton distributions for various colliding systems with 
double power law formula: 

Sa(dN ch /d V , v/iii) = a x (dN ch /d V f x (y^E) 7 (13) 

The fitting results are shown in the Table. [TJ The Sa can be well described by the power law 



Table 1. Fitting parameters for Sa as a function of dN c h/dn 



Parameters 


Value 


Approx. 


X z /ndf 


34.08/37 


0.92 


a 


3.669+0.1358 


11 


/3 


0.0853+0.003927 


1 
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7 


-0.6796+0.00693 





formula Sa = — x dif* ) 12 > wnere the s is the square of the center of mass energy. For high 
energy heavy ion collisions the temperature is approximately constant and the ratio ^b/T « 1, 



thus we have the approximation /Ub/T ~ tanh(^g/r) = Sa = ^ x (-^ d ^ h ) 12 . This denotes 
the relation between hb/T and the charged particle density and colliding energy for high energy 
nuclear collisions. 

Multiplicity fluctuations and inclusive yields are two basic properties in high energy heavy 
ion collisions. For a thermal system, both the fluctuations and yields should be described by 
the thermodynamic parameters (hb and T), which completely determine the properties of the 
thermal system. The fluctuation observable Sa of most central net-proton distributions versus 
thermodynamic parameter hb/T ratios, which are extracted from the thermal model fit of the 
particle ratio, is shown in the Fig. [8] for various colliding systems. Lattice QCD calculations 
with N T = 6 and T c = 175 MeV and the HRG model relation Sa = tanh(/iB/T) are also 
plotted in Fig. [8] for comparison. We find that high energy heavy ion collisions, such as Au+Au 
and Cu+Cu collisions, are consistent with Lattice QCD and HRG model calculations, while 
the elementary p+p collision deviate from the HRG model calculations. In addition to the 
perfect description of the particle yields by the thermal model, the agreement of higher order 
fluctuations with thermal model predications provide further evidence that the colliding system 
has achieved thermalization in most central high energy heavy ion collision. 

5. Summary 

The higher moments of net-proton multiplicity distributions measured in heavy ion collision 
experiment have been applied to search for the QCD critical point, due to the high sensitivity 
to the correlation length (£). The beam energy and system size dependence for higher moments 
(M, a, S, k) as well as moment products (na 2 , Sa) of net-proton multiplicity distributions have 
been presented with a broad energy range and different system sizes, which include Au+Au 
collisions at ^/s^ = 200, 130, 62.4, 39 and 19.6 GeV, Cu+Cu collisions at y /s^ = 200, 62.4 and 
22.4 GeV, d+Au collisions at ^/s^ = 200 GeV, p+p collisions at ^/s^ = 200 and 62.4 GeV. 
The moment product Ka 2 shows no centrality dependence while Sa shows a weak centrality 
dependence. The energy dependence is studied by comparing the results from the Au+Au 
200 GeV to those from the BES energies. The moment products na 2 and Sa of net-proton 
distributions from most central Au+Au collisions are consistent with Lattice QCD and HRG 
model calculations at high energy (200, 130, 62.4 GeV) while deviating from (smaller than) 
HRG model calculations at ^/s NN = 39 GeV. Lattice QCD calculations show negative value for 
Ka 2 at ^/s NN = 19.6 GeV. The deviations could potentially be linked to chiral phase transitions 
and QCD critical point. But the experimental data is with large error bar due to the limited 
statistics. Fortunately, this ambiguity can be clarified soon by 19.6 GeV data taken in year 2011 
with higher statistics. Recent model calculations show that the na 2 value will always be smaller 
than its Poisson statistical value 1, when QCD critical point is approached from the high energy 
cross-over side. 

On the other hand, the mutual agreements between the /j-b/T extracted from thermal model 
fits of particle ratio and from the event-by-event fluctuations observable Sa of net-proton 
distributions provides further evidence of thermalization of the hot dense matter created in 
the heavy ion collisions. Further, the na 2 and Sa of net-proton distributions are consistent with 
Lattice QCD and HRG model calculations at high energy, which support the thermalization 
of the colliding system. The deviations of the Ka 2 and Sa of net-proton distributions for 
Au+Au central collisions from HRG model predications at low energies are not well understood. 
This may result from the non-applicability of grand canonical ensemble or the appearance of 
QCD critical point and chiral phase transitions at low energies. However, it should be further 
investigated. 
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